rm(list=ls())
setwd("/Volumes/Disk2/Future_PM/Dataverse")
library(fields); library(maps); library(ncdf4);library(abind)
source('Function/get_geo.R')
source('Function/get_met.R')
source('Function/functions.R')

load("Figures/Figure 3/monthly_CV.Rdata")
mai=c(0.1, 0.1, 0.2, 0.1)
mgp=c(1.2, 0.4, 0)
tcl=-0.2
ps=12

close.screen(all.screens = TRUE)
m <- rbind(c(0, 0.5, 0.5, 0.9), c(0.5, 1, 0.5, 0.9),c(0.2,0.8,0.44,0.7))
split.screen(m)

cr=0.4
screen(1)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
ap=find.cor(monthly_PM.std, local.PM,plim=1)
image(lon.frame,lat.frame, ap^2,zlim=c(0,0.70),col=tim.colors(32),xaxt='n',yaxt='n',xlab='',ylab='');
map("world", add = T)
clines=contourLines(lon.frame, lat.frame, ap^2,level=c(cr))
for(k in 1:length(clines)){
	# lines(clines[[k]]$x,clines[[k]]$y,lwd=1,col=1,lty=2)
}
title('(a) Local meteorology',cex.main=0.9,font.main=1)
text(x=-70,y=26,c('34%'),font=1,cex=1.0)

screen(2)
par(mai=mai, mgp=mgp, tcl=tcl, ps=ps)
ap=find.cor(monthly_PM.std, hybrid.PM,plim=1)
image(lon.frame,lat.frame, ap^2,zlim=c(0,0.70),col=tim.colors(32),xaxt='n',yaxt='n',xlab='',ylab='');
map("world", add = T)
clines=contourLines(lon.frame, lat.frame, ap^2,level=c(cr))
for(k in 1:length(clines)){
	# lines(clines[[k]]$x,clines[[k]]$y,lwd=1,col=1,lty=2)
}
title('(b) Local meteorology + synoptic factors',cex.main=0.9,font.main=1)
text(x=-70,y=26,c('43%'),font=1,cex=1.0)

screen(3)
mai=c(0.1, 0.1, 0.05, 0.2)
par(mai=mai, mgp=mgp, tcl=tcl, ps=9)
image.plot(ap^2,legend.shrink=0.9,legend.only=TRUE,zlim=c(0,0.70),col=tim.colors(32),horizontal=TRUE,legend.width=2.5,axis.args = list(mgp = c(0, 0.5, 0),tcl=-0.2,padj=-1
))
text(expression(paste('R'^'2','',sep='')),x=par("usr")[2]-0.23,y=par("usr")[3]+0.51,srt=0, adj = 0,xpd=TRUE,cex=1.2)
